#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_sf_legendre.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>
#include <iostream>
#include <string>
#include <time.h>     
     int main (){
  time_t start_time,end_time;
  start_time = time(NULL);//Record time at start of program.
static const double pi = 3.141592653589793238462643383279502884197; // Archimedes constant pi

  //Convert coordinates to radians.
  double lat1 = 90.0*pi/180;
  double lon1 = 0.0*pi/180;
  double lat2 = 72.0*pi/180;
  double lon2 = 0.0*pi/180;
  
  //Calculate central_angle (in radians) from Haversine formula.
  double central_angle = 2*asin(sqrt(pow(sin((lat1-lat2)/2),2)+cos(lat1)*cos(lat2)*pow(sin((lon1-lon2)/2),2)));
  
  printf("distance(km): %f",central_angle*6372.795477598);



  //printf("time: %ld",(long)start_time);

/*
            int gsl_fft_halfcomplex_unpack (const double halfcomplex_coefficient[], gsl_complex_packed_array complex_coefficient, size_t stride, size_t n)

    This function converts halfcomplex_coefficient, an array of half-complex coefficients as returned by gsl_fft_real_transform, into an ordinary complex array, complex_coefficient. It fills in the complex array using the symmetry z_k = z_{N-k}^* to reconstruct the redundant elements. The algorithm for the conversion is,


            
              //To perform an FFT on a vector argument, such as gsl_vector_complex * v, use the following definitions (or their equivalents) when calling the functions described in this chapter:

     //gsl_complex_packed_array data = v->data;
     //size_t stride = v->stride;
     //size_t n = v->size;*/



       int i;//, n = 20;//change below too!  why was 'n' okay in gslref?
       double dummy;
       long num=50000;
       double data[50000];

  double x=0.5;
  
  //choice = 0, one at a time
  //choice = 1, arrays
  int choice = 1;
  if(choice==0) for(i=0;i<num;i++) dummy = gsl_sf_legendre_Pl(i, x);//printf("P%d(%f) = %f\n",i,x,gsl_sf_legendre_Pl (i, x));
  else if (choice == 1){
    gsl_sf_legendre_Pl_array(num, x, data);
  }
  
  end_time = time(NULL);//Record time at end of program.
  printf("Program took %d seconds for %ld Legendre calcs\n",(int)(end_time-start_time),num);

     /*
       gsl_fft_real_wavetable * real;
       gsl_fft_halfcomplex_wavetable * hc;
       gsl_fft_real_workspace * work;
       
       //I've got to initialize this otherwise I get a warning- no other reason
       gsl_complex_packed_array complex_coefficient=data;
       
       for (i = 0; i < n; i++)
         {
           data[i] = 0.0;
         }
     
       for (i = n / 3; i < 2 * n / 3; i++)
         {
           data[i] = 1.0;
         }
     
       for (i = 0; i < n; i++)
         {
           //printf ("%d: %e\n", i, data[i]);
         }
       //printf ("\n");
     
       work = gsl_fft_real_workspace_alloc (n);
       real = gsl_fft_real_wavetable_alloc (n);
     
       gsl_fft_real_transform (data, 1, n, 
                               real, work);

       gsl_fft_halfcomplex_unpack (data, complex_coefficient, 1, n);
       //printf ("Original FT, before lowpass filter: \n");
       for (i = 0; i < n; i++)
         {
           //printf ("r[%d]: %e, i: %e, sqrtsumsqr: %e\n", i, complex_coefficient[2*i], complex_coefficient[2*i+1], sqrt(pow(complex_coefficient[2*i],2)+ pow(complex_coefficient[2*i+1],2)));
         }
       //printf ("\n\n");
     
       for (i = 11; i < n; i++)
         {
           data[i] = 0;
         }
     
       hc = gsl_fft_halfcomplex_wavetable_alloc (n);
     
       gsl_fft_halfcomplex_inverse (data, 1, n, 
                                    hc, work);

       gsl_fft_real_transform (data, 1, n, 
                               real, work);

       gsl_fft_halfcomplex_unpack (data, complex_coefficient, 1, n);
       //printf ("FT after lowpass filter: \n");
       for (i = 0; i < n/2; i++)
         {
           //printf ("r[%d]: %e, i: %e, sqrtsumsqr: %e\n", i, complex_coefficient[2*i], complex_coefficient[2*i+1], sqrt(pow(complex_coefficient[2*i],2)+ pow(complex_coefficient[2*i+1],2)));
         }
       //printf ("\n\n");

     
       gsl_fft_real_wavetable_free (real);

       gsl_fft_halfcomplex_wavetable_free (hc);
     
       for (i = 0; i < n; i++)
         {
           //printf ("%d: %e\n", i, data[i]);
         }
     
       gsl_fft_real_workspace_free (work);
       
       //printf("%20s%20s\n","X value","Y value");
       //printf("%20.8e%20.8e\n",1.666666666666,2.333333333333);//*/
       return 0;
     }
